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Use of edge-based finite elements for 

SOLVING THREE DIMENSIONAL SCATTERING 

PROBLEMS 


A. Chatterjee, J.M. Jin and J.L. Volakis 
Radiation Laboratory 

Department of Electrical Engineering and Computer Science 

University of Michigan 
Ann Arbor, MI 48109-2122 

Abstract 

Edge-based finite elements are free from the drawbacks associated with node- 
based vectorial finite elements and are, therefore, ideal for solving three- 
dimensional scattering problems. The finite element discretization using edge 
elements is checked by solving for the resonant frequencies of a closed in- 
homogeneously filled metallic cavity. Significant improvements in accuracy 
are observed when compared to the classical node-based approach with no 
penalty in terms of computational time and with the expected absence of 
‘spurious’ modes. A performance comparison between the edge-based tetra- 
hedra and rectangular brick elements is carried out and tetrahedral elements 
are found to be more accurate than rectangular bricks for a given storage 
intensity. A detailed formulation for the scattering problem with various 
approaches for terminating the finite element mesh is also presented. 



1 Introduction 


The problem of computing the electromagnetic field scattered by a material 
body, when a field is incident on it, has great theoretical and practical im- 
portance. It is, usually, desirable to relate the characteristics of the scattered 
field in the far zone with the shape of the scatterer. Finite elements can be 
used to model complicated geometries but are of no help in computing the 
far field. The far field effect is thus represented in finite element analysis 
by enclosing the target within a fictitious boundary and enforcing a special 
boundary condition, often called the absorbing boundary condition, on the 
outer boundary. In this report, we discuss the choice of finite elements in the 
inner region. As a first step, we solve Maxwell’s equations for the resonances 
of a closed cavity to verify the finite element discretization and decide on the 
type of basis functions and elements to use. 

The occurrence of ‘spurious’ modes [1] in the node-based finite element 
approach has plagued the computation of cavity eigenvalues for years. This 
difficulty can be circumvented with the introduction of a penalty term[2] to 
render the finite element vector field solutions non-divergent. However, it 
is difficult to satisfy continuity requirements across material interfaces[3] or 
for geometries with sharp edges[4] using classical finite-elements, obtained 
by interpolating the nodal values of the vector field components. This is es- 
pecially critical in solving scattering problems using classical finite elements 
where special techniques need to be used for satisfying normal and tangential 
continuity of the electric displacement and the electric field, respectively [5]. 
Edge elements, a type of vector finite element with their degrees of free- 
dom associated with the edges of the mesh, have been shown to be free of 
these shortcomings^], [7]. Edge-based finite elements are, therefore, a nat- 
ural choice for treating three dimensional geometries. Generally these lead 
to more unknowns but the higher variable count is balanced by the greater 
sparsity of the finite element matrix so that the computation time required 
to solve such a system iteratively with a given accuracy is less than the 
traditional approach[8|. 

In this report, we have discussed the formulation for both the cavity and 
the scattering problems but have presented results only for the former. We 
have solved for the eigenvalues of an arbitrarily shaped metallic cavity using 
node-based and edge-based vector finite elements. The computed data are 
then compared with analytical results for empty and partially filled cavities. 
A comparison between the storage intensity and computational accuracy for 
edge-based rectangular bricks and tetrahedra is also presented. Finally, we 
compute the eigenvalues of a metallic cavity with a ridge along one of its 
faces. 
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2 Formulation 


2.1 Derivation of finite element equations 

Consider a three dimensional inhomogeneous body occupying the volume V. 
To discretize the electric field E inside this volume, we subdivide it into a 
number of small tetrahedra or rectangular bricks, each occupying the volume 
V e (e = 1, 2, • • • , M), where M is the total number of elements. Within each 
element, the electric field satisfies the vector wave equation 

Vx-VxE-^ r E = 0 (1) 

Hr 

in which fi T is the permeability of the medium, t T is the medium permittivity 
and k 0 is the free space wave number. For a numerical solution of ( 1), we 
expand the electric field within the eth volume element as 

m 

E = £ £'W* (2) 

j=i 

where W' are edge-based vector basis functions, E \ ■ denote the expansion 
coefficients of the basis, m represents the number of edges in the element and 
the superscript stands for the element number. On substituting ( 2) into ( 1) 
and upon applying Galerkin’s technique, we obtain 

E E j J y W i • ( V X y V X W i “ ^ £ r W j) dv = 0 


Further, upon making use of basic vector identities and the divergence the- 
orem, we obtain the weak form of Maxwell’s equation 



— (V x W?) . (V x Wj) - k 2 0 e r W'.W' 


dv = jk 0 Z 0 W'.(n x H )ds 
J S c 

(3) 


where n x H is the tangential magnetic field on the exterior dielectric sur- 
face, S e denotes the surface enclosing V e and Z 0 is the free-space intrinsic 
impedance. Equation ( 3) can be conveniently written in matrix form as 

mpn = kl\B'}{E’} + {C*) (4) 

where 

A 'j = I — (V X WJ) . (V x W') 

JV e fi r 

B\ = / e r W*.W> 

JV' 3 

C? = jk o Z 0 l W*.(n x H )ds 

J S e 

On assembling the equations from all the elements making up the geometry, 
we obtain the system 


( 5 ) 

( 6 ) 
(7) 
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( 8 ) 


M A/ M 

e=l e=l e=l 

Due to the continuity of tangential H at the interface between two di- 
electrics, an element face lying inside the body does not contribute to the last 
term of ( 8) in the final assembly of the element equations. As a result, the 
last term of ( 8) reduces to a column vector containing the surface integral 
of the tangential magnetic field only over the exterior surface of the body. 
In this paper, we are interested in a perfectly conducting surface surround- 
ing the volume V and in this case, the surface integral vanishes altogether. 
Therefore, this system can be more compactly written as 

|A]{£} = A [£]{£} (9) 

where [A] and [B] are N x N symmetric, sparse matrices with N being the 
total number of edges resulting from the subdivision of the body, {£} is a 
N x 1 column vector denoting the edge fields and A = k% gives the eigenvalues 
of the system. Solving ( 9), we obtain the resonant wavenumber k 0 and the 
resonant field distribution {E} as well. 

In the case of the scattering problem, ( 3) can be conveniently written in 
matrix form as 

A' e ]{E'} = {C e } 


where 


A' e 

A ij 

= i (V x Wf) . (V x W‘) - fcJe.Wf.w;) dv 

•'■'e \f*r / 

(11) 

c t 

= ]k 0 Zo f Wf.(n x H )ds 

J Sc 

(12) 

Summing over all M elements in the geometry, we obtain a system of equa- 
tions whose solution yields the field components over the entire body. 


M M 

EK]{£'} = £{<?*} 

t-\ e=l 


which gives 

[A'J{£} = (C) 

(13) 


where {i?} are defined as in ( 9), [A'] is a N x N matrix and {C} is a 
column vector of size N related to the tangential magnetic field over the 
exterior surface of the body. On separating the on-surface edges from the 
edges inside the body for ease of representation, we obtain 

[A'„H£.}+[A'„){£,} = {C,} 

[Ay {£,} + [/!'„]{£,} = 0 

{C.} = [£'„]{£,} (14) 


( 10 ) 
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where the subscript s denotes the edges on the surface and i represents the 
edges inside the body. It is thus readily seen that ( 14) relates the electric field 
inside and on the surface of the body to the on-surface tangential magnetic 
field. However, ( 14) contains two unknown vectors { E } and { H } and the 
system can, therefore, be solved only when we relate the electric field with the 
magnetic field. This is done by introducing an absorbing boundary condition 
at the outer boundary of the mesh as discussed in section 2.3. 


2.2 Basis functions 


The vector edge basis functions for rectangular bricks are outlined in detail 
in [9]. The basis functions for tetrahedral elements merit detailed mention 
since they are more involved. 

Vector fields within tetrahedral domains in three dimensional space can be 
conveniently represented by expansion functions that are linear in the spatial 
variables and have either zero divergence or zero curl. The basis functions 
defined in [8] are associated with the six edges of the tetrahedron and have 
zero divergence and constant curl. Assuming the four nodes and the six 
edges of a tetrahedron are numbered according to Table 1, the vector basis 
functions associated with the (7 — i)th edge of the tetrahedron are defined as 


with 


W 7 _, = 


f 7 _, + g 7 _, x r, r in the tetrahedron 
0, otherwise 


h-i 

g7 — . 


h-i 

W r " x r ” 

b i^7—i Cj 

6V 


(15) 

(16) 
(17) 


where i = 1, 2, • • • , 6, V is the volume of the tetrahedral element, e, = (r,- 2 — 
r.J/6, is the unit vector of the ith edge and 6,- = jr tJ — r M | is the length of 
the ith edge. All distances are measured with respect to the origin. 

Since there are two numbering systems, local and global, a unique global 
edge direction is defined (e.g., always pointing from the smaller node number 
to the larger node number)[10] to ensure the continuity of n x E across all 
edges. This implies that ( 15) should be multiplied by (-1) if the local edge 
vector (as defined in Table 1) does not have the same direction as the global 
edge direction, i.e. the local vector points from the larger node number to 
the smaller node number. Finally, since V- W ,• = 0 the electric field obtained 
through ( 2) exactly satisfies the divergence equation within the element, i.e. 
V • E = 0. Therefore, the finite element solution is free from contamination 
of spurious solutions. 


2.3 Mesh termination 

Differential equation methods, such as finite elements, can only solve bound- 
ary value problems. Since electromagnetic problems are open boundary- 
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infinite domain types, a means to truncate the solution domain to lie within 
a finite boundary must be found. On this boundary, a condition is enforced 
thus ensuring that the fields will obey the Sommerfeld radiation condition 
at distances asymptotically far from the object. These absorbing boundary 
conditions (ABCs) have a significant advantage over the global methods of 
solving unbounded problems using finite elements in that they are local in 
nature. Due to this, the sparse matrix structure of the finite element formu- 
lation is retained. One disadvantage, however, is that ABCs are approximate 
and do not model the exterior field exactly. 

The objective of absorbing boundary conditions is to truncate the finite 
element mesh with boundary conditions that cause minimum reflections of 
an outgoing wave. These ABCs should provide small, acceptable errors while 
minimising the distance from the object of interest to the outer boundary. 
This minimal distance is required to reduce the number of unknowns in the 
problem for computational efficiency. 

There are two ways to approach this problem. The first approach is to 
employ a three dimensional vector boundary condition for terminating the 
finite element mesh of the body described in section 2.1. The second method 
involves placing an artificial conducting boundary (electric or magnetic) to 
truncate the infinite region and then coating the inner surface of the bound- 
ary with a layer or several layers of fictitious dielectric whose thickness and 
constitutive parameters are chosen to minimise the non-physical reflections 
of the scattered field over a wide range of incidence angles[ll]. 


2.3.1 Vector ABC 


We begin with the Wilcox representation [12] of the electric field which has 
an expansion 


E( r ) = e Jfcr ^ 


r " r n 

n= 0 


From ( 18), we get 


V xE={„x^}e4^ 

n=l 


(18) 


(19) 


where A nt = r x A n is the transverse component of A n and, for a vector F, 
jDjF is given by 

1 


D x F = 


sinO 

1 


'!<*"*> - % 


- sinOF * 

0 + 

\f> - a ;;i 

d0 


de 


sin9 

Using the recursion relation 

-2jknA nt = n(n - l)A n _ 1;( + D A A n _ x 


4 > 


(20) 
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where 


D 4 A n 


Dg A n 


D 4 , A n 


(D A e n + D e A n ) + ^ A n ) 

5A; 1 fi 2cosfldA* 

sin 2 6 n sin 2 6 dO 

J_dK _ _L_a* 4- 2c05ga/L " 

smfl sin 2 6 n sin 2 9 d(f> 


and D is Beltrami’s operator[10], we can derive the representation correct to 
r -4 . Applying the recursion relation in ( 19) yields the desired relationship 
for the vector ABC: 


V x E = a(r)E + P(r)D 4 E 


( 21 ) 


where 


a(r) 

P( r ) 




_1 1 

2 jkr 2 (1 + l/jkr) 



( 22 ) 

(23) 


The ABC formulated above was derived for spherical boundaries and hence 
would be storage intensive and numerically inefficient when used to terminate 
the mesh of long and thin geometries. It would be highly desirable to choose 
an outer boundary that conforms to the shape of the object. It is the authors’ 
intention to extend the boundary condition discussed above to conformal 
boundaries. Another way to introduce a conformal outer boundary is to use 
the fictitious absorber model discussed below. 


2.3.2 Fictitious absorber model 

In this section, our intention is to develop a model which can absorb the 
reflections of the scattered field from the artificial outer boundary over as 
wide a range of incidence angles as possible. The reflection coefficient on 
applying the appropriate boundary condition can be symbolically written as 

A y(Ui^rl>/^rlj^2 5^r2)/^r27'**j^n>^rn>^rni^}^) (24) 

for an n-layer absorber. Here, and /i r , denote the thickness, relative 

permittivity and relative permeability of the ith dielectric layer. Using a 
multi-dimensional minimization algorithm such as that based on the downhill 
simplex method, we can find a set of t,, e rt , /i r , that minimises the magnitude 
of the reflection coefficient over a wide range of incidence angles. Further, 
the outer boundary of the mesh can be made conformal to the shape of the 
scatterer reducing the number of unknowns compared to the earlier approach 
which works solely for spherical outer boundaries. 
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2.4 Solution of the finite element equations 

2.4.1 Vector ABC 

An inspection of ( 13) reveals that for an inhomogeneous body, there is no a 
priori information about the tangential magnetic field over the exterior sur- 
face of the body. Relation ( 13) therefore contains two unknown vectors, { E } 
and {C}, and thus another condition is required involving the two variables 
to permit an evaluation of the fields inside and on the surface of the body. 
This condition relating the tangential electric field to the tangential magnetic 
field on the surface is provided by ( 21). Since the ABC in ( 21) refers to the 
scattered field, we can rewrite it as 

VxE; = a(r)EJ + P(r)D<E° 

= -i- [<*(r)E* + /?(r)ft,E*,] 

to fl 

= (25) 

where K. = ^ [<*(r) + (3(r)D 4 \ and the subscript s denotes the field on the 
surface and the superscript s represents the scattered field. As the total field 
is a sum of the incident field and the scattered field, therefore from ( 25), we 
obtain 


H,-H‘ nc = /C(E S -E’ nc ) (26) 

Since {# s } = [i? JS ] -1 {C' J } from ( 14), we obtain on substituting ( 26) into 

(14) 

{[A'„]-[B„]K:){E.} + [A' ri ){E i } 

[A'i.] {E s } + [A'u] {E,} 

The above equation can thus be solved 
inside and on the surface of the body. 

2.4.2 Fictitious absorber model 

Since an artificial outer boundary is placed in this case, the column vector 
{C} in ( 13) is a function of the incident field and the integration is carried 
out over the outer and inner surfaces. Therefore, we can rewrite ( 14) as 

K..1 (£» + Kd {£?} = {c"“} 

K'..l (£,'} + l-4".w] {£*} = 0 (28) 

where the superscript inc stands for the incident field, the superscript s 
denotes the scattered field, the subscripts are as defined in ( 14) and [ A"] 
contains contributions from the surface integral as well as \A'l. 


= [bj ({#;■“} - k {£r}) 

= 0 (27) 

for the unknown electric fields both 
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3 Results 


In Table 2, the eigenvalues for a lcmx .5cmx.75cm rectangular cavity are 
calculated using nodal basis functions. This particular derivation is outlined 
in [13] and the existence of ‘spurious’ modes is clearly observed when the 
zero divergence condition is not rigidly enforced (s = 0). For non-zero values 
of s, ‘spurious’ modes do not appear. However, arbitrary values of s lead to 
inaccuracies and, therefore, the most suitable value needs to be determined 
by trial and error. 

Using edge- based rectangular bricks results in better accuracy as observed 
in Table 3. However, on comparing bricks with tetrahedra, it is seen from 
Table 4 that the tetrahedral edge element is the best in terms of accuracy for 
the same storage requirement. The maximum edge length for the rectangular 
brick elements was .15cm whereas that for the tetrahedral elements was .2cm. 

In Tables 4-7, we compare the eigenvalues computed using edge-based 
tetrahedral finite elements with their analytical counterparts. The finite ele- 
ment mesh was generated using SDRC I-DEAS, a commercial pre-processing 
package and it is seen that the numerical results are in excellent agreement 
with the exact values for both homogeneous and inhomogeneous cavities. In 
Table 4, the edge-based approach predicts the first six distinct non- trivial 
eigenvalues with less than 4 per cent error for an empty rectangular cavity. 
The exact eigenvalues of the half-filled cavity as described in Table 5 are 
computed by solving the transcendental equation obtained by matching the 
tangential electric and magnetic fields at the air-dielectric interface. As seen, 
these results agree with the results of the finite element code to within 1 
per cent. Table 6 compares the analytical and computed eigenvalues of a 
cylindrical cavity having a base radius of 0.5cm and a height of 0.5cm. The 
eigenvalues are again seen to be quite accurate and the accuracy is expected 
to increase further with higher mesh density. Similar comparisons are given 
in Table 7 of a sphere of Icm radius. Finally, Table 8 presents the eigenvalues 
of the geometry shown in Fig. 1. This is basically a closed metallic cavity 
with a ridge along one of its faces. We note again that it is difficult to treat 
such a geometry having sharp edges and corners using node-based finite el- 
ements. The storage requirement for the geometries discussed above can be 
further reduced by taking symmetry into account. Since we are interested 
in developing a code for arbitrary geometries, boundary conditions taking 
advantage of symmetry have not been introduced. 

It is noted that as the degeneracy of the eigenvalues increases, the matrix 
becomes increasingly ill-conditioned and the numerical solution correspond- 
ingly less accurate[14]. This is clearly observed from the data in Table 7 for 
the case of a perfectly conducting hollow spherical cavity. Since the second 
lowest TM mode has five-fold degeneracy, the computational error is seen 
to be the greatest. However, for the partially filled rectangular cavity, the 
absence of degenerate modes gives results which are accurate to within 1 per 
cent of the exact eigensolutions. 


8 



4 Conclusions and future work 


It has been shown that the resonant frequencies of an arbitrarily shaped 
inhomogeneously filled metallic resonator can be computed very accurately 
using edge-based tetrahedral elements. The nodal method is not very reliable 
and is dependent on the value of the s parameter which is a measure of how 
strongly one wants to enforce the divergence-free condition. Edge-based rect- 
angular bricks do not provide as good an accuracy as edge-based tetrahedral 
elements and is limited to a small class of geometries. It has, therefore, been 
established that it is best to use edge-based tetrahedral elements for solving 
three-dimensional problems. 

In the future, we intend to work on three dimensional scattering problems 
using edge-based tetrahedra and employing the mesh termination techniques 
discussed in section 2.3. We also plan to focus on the derivation of new 
absorbing boundary conditions conformal to the shape of the scatterer. 


5 Appendix 

5.1 Derivation of matrix elements 

The derivation of the matrix elements in ( 9) amounts to evaluating the 
integrals in ( 5) and (6). Therefore, from ( 5), we have 

/-(VxW').(VxWp = -g„ gi K (29) 

JV C fi T J fi r 

since V x W* = 2g t . The evaluation of the integral in ( 6) is more cumber- 
some. Substituting into ( 6) the basis functions defined in ( 16) and ( 17), 
we obtain 

«, / W'.W )iv = e, / {(fi-fj) + (r.D) + (g, x r).(g, x r)} dv (30) 

•JV C JV C 

= £r (A + I2 + I3) 

where 

D = ( f i x g j) + ft x g.) 

and 

h = / v f,f idv (31) 

h = J y r.D dv (32) 

13 = Jv ^ S ' X r ^ Sj X T ^ dv ^ 
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Since f is a constant vector, /j reduces to 


h = fi-fi V e (34) 

Since 

4 

x = ) ] L{X{ 

«=i 

4 

y = ]£ 

t=i 

4 

z - 51 

t=l 


where L, are the shape functions for the tetrahedral element and x,, y,, z,(i = 
1, • - • ,4) denote the x,y and z co-ordinates of the vertices of the tetrahedral 
element. Using the standard formula for volume integration within a tetra- 
hedral element and simplifying, we have 


h 


k 

4 


4 4 4 

D x Y Xi : + D y Y Vi + Dz Y z i 

1=1 i=i t=i . 


(35) 


where D m is the mth component of D. The evaluation of / 3 can be simplified 
by the use of basic vector identities. Therefore, 


h = 


g. gj J v M 2 dv ~ J v (g«- r ) (gj- r ) dv 

{diy9jy T 9>z9jz) I x dv -j- ( 9ix9jx T 9%z9jz) I V dv -f- {9ix9jx T 9iy9jy) I z dv 
JVt JV* ^ Ve 

( 9ix9jy T 9jx9iy) I xydv — {9ix9jz 9jx9iz ) / zxdv ( 9iy9jz T 9jy9iz) I yzdv 
J V e J V c J v c 

(36) 


where gi m represents the mth component of the vector g*. Each of the inte- 
grals in ( 36) can be easily evaluated analytically and the result expressed in 
the following general form: 



aia m dv 


Ve 

20 


' 4 
.1=1 


4 4 

mi y ] Q-li y a mi 

t=l i=l 


(37) 


where /, m = 1 , - - - , 3 and ai represents the variable x, a 2 stands for the 
variable y and a 3 denotes the variable z. 
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Edge no. 


*2 

■ 

i 

2 

fKwj 

i 

3 

■ 

i 

4 


2 

3 


4 

2 


3 

4 


Table 1: Tetrahedron edge definition 


Mode no. 

Exact 

s=0 

s=.5 

s=l 

s=5 

spurious 


2.010 




spurious 


3.312 




spurious 


4.443 




spurious 


4.928 




TE101 

5.236 

5.318 

5.427 

5.447 

5.571 

spurious 


5.759 




spurious 


6.255 




spurious 


6.804 




TM 110 

7.025 

7.067 

6.574 

7.543 

7.940 

spurious 


7.305 




TE Oil 

■jgSjj 

7.389 

7.476 

7.635 

8.014 

TE 201 


7.640 

7.659 

8.072 

8.348 


Table 2: Eigenvalues for an empty 1cm x0.5cm x0.75cm rectangular cavity 
computed using nodal basis (230 unknowns) 


Mode no. 

Analytical 

Computed 

Error(%) 

TE 101 

5.236 

5.307 

-1.36 

TM 110 

7.025 

7.182 

-2.23 

TE Oil 


7.725 

-2.58 

TE 201 

| . 

7.767 

-3.13 

TM 111 

8.1789 

8.350 

-2.09 

TE 111 

8.1789 

8.350 

-2.09 

TM 210 

8.886 

9.151 

-2.98 

TE 102 

8.947 

9.428 

-5.38 


Table 3: Eigenvalues for an empty 1cm x0.5cm x 0.75cm rectangular cavity 
using edge-based rectangular bricks (270 unknowns) 
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Mode no. 

Analytical 

Computed 

Error(%) 

TE 101 

5.236 

5.213 

.44 

TM 110 

7.025 

6.977 

.70 

TE Oil 



1.00 

TE 201 


. c ! 3 

-.56 

TM 111 

8.1789 

7.991 

2.29 

TE 111 

8.1789 

8.122 

.70 

TM 210 

8.886 

8.572 

3.53 

TE 102 

8.947 

8.795 

1.70 


Table 4: Eigenvalues for an empty 1cm x0.5cm x0.75cm rectangular cavity 
(260 unknowns) 


Mode no. 

Analytical 

Computed 

Error(%) 

TEz 101 

3.538 

3.534 

.11 

TEz 201 

5.445 

5.440 

.10 

TEz 102 

5.935 

5.916 

.32 

TEz 301 

7.503 

7.501 

.04 

TEz 202 

7.633 

7.560 

.97 

TEz 103 

8.096 

8.056 

.50 


Table 5: Eigenvalues for a half-filled 1cm xO.lcm xlcm rectangular cavity 
(e r = 2) (192 unknowns) 


Mode no. 

Analytical 

Computed 

Error(%) 

TM 010 

4.810 

4.809 

.02 

TE 111 

7.283 


1.1 



7.288 

-.07 

TM 110 

7.650 


.22 




-.97 

TM 011 

7.840 

7.940 

-1.28 

TE 211 

8.658 

8.697 

-.45 



8.865 

-2.39 


Table 6: Eigenvalues for an empty cylindrical cavity of base radius 0.5cm 
and height 0.5 cm (380 unknowns) 
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Mode no. 

Analytical 

Computed 

Error(%) 

TM 010 

2.744 

2.799 

-2.04 

TM 111, even 


2.802 

-2.11 

TM 111, odd 


2.811 

-2.44 

TM 021 

3.870 

3.948 

-2.02 

TM 121, even 


3.986 

-2.99 

TM 121, odd 


3.994 

-3.20 

TM 221, even 


4.038 

-4.34 

TM 221, odd 


4.048 

-4.59 

TE Oil 

4.493 

4.433 

1.33 

TE 111, even 


4.472 

.47 

TE 111, odd 


4.549 

-1.25 





























The following is a listing of the computer programs 
generating the node-based and edge-based elements, 
interfacing with SDRC IDEAS and computing the eigenualues of an 
inhomogeneously filled, arbitrarily shaped cauity. 


i 
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Processes a universal file obtained from IDEAS and converts the nodal info 
to edge info needed for constructing an edge-based three-dimensional finite 
element mesh using tetrahedral elements. 

Stores the node numbers and respective nodal coordinates in 'noddy' 

Stores the element numbers and corresponding nodes in 'elno' 

Processes 'noddy' and 'elno' and stores the edge nos. and nodal connections 
in 'glob' and edge nos. with corresponding edge vector in 'edgy' 

Note: Edge overlaps are taken care of. 

Storage limit: 800 nodes, 3000 elements, 4300 edges 


program unv_f ile__processor 
character string*80, yastrn*40, unv*20 

integer nl (3000,4) >tab(4300,2) ,nn(100) ,tr(3000) , nun (3000) , 
1 be (3000, 3) ,edgv (18000, 3) , mat (3000) 

real x (800) ,y (800) , z (800) 
common /bank/nl, x, y , z, ne, mat 
common /dbase/edgv 
common /local/ncount 

write ( 6, * ) ' 1 ) Cavity problem 2) Scattering problem' 
read(5, *) ivar 

write (6, *) 'Name of universal file ?' 

read (5, ' (A) ' ) unv 

open (1, file=unv) 

open (2, f ile=' enoddy' ) 

open (3, f ile=' elno' ) 

open (4, f ile=' eglob' ) 

open (7, file=' edgy' ) 

open (8, f ile=' esurfed' ) 

.Processing universal file for info on various parameters 
read(l, ' (A) ' ) string 
1=0 

if (string (4 : 6) . eq. ' 151' ) then 
write (6, *)' Encountered header' 
go to 15 

elseif (string (4 : 6) . eq. ' 747' ) then 
read(l,' (A) ') string 
if (string (5 : 6) . ne . ' -1' ) then 
if (string (9 : 9) . ne . ' ') then 
write (4, *) string (2 : 13) 
do ik=l , 4 

readd, ' (A) ' ) string 
enddo 
endif 
go to 2 
endif 

.Processing nodal info; data stored in 'noddy' 
elseif (string (4 : 6) . eq. ' 15') then 
read(l,' (A)') string 
if (string (5 : 6) . ne . ' -1 ' ) then 

write (2, *) string (7 : 10) , ' ', string (41 : 53) , ' ', 

1 , string (67 : 79) 

1=1+1 
go to 5 
endif 
11=1 

write (6, *) 'There are ',1,' nodes.' 
write (4, *) 1 

.Processing element info; data stored in 'elno' 


string (54 : 66) 


t t 
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elseif (string (4 : 6) . eq. ' 71') then 
read(l, ' (A) ' ) string 
if (string (5 : 6 ) . ne . ' -1 ' ) then 
read(l,' (A)')yastrn 

write (3, *) string (7 : 10 ) , ' ' , yastrn (7 : 10) , ' ' , yastrn ( 17 : 20 ) , 

1 ' ' ,yastrn (27 : 30) , ' ', yastrn (37 : 40) , ' ', string (49 : 50) 

1=1+1 
go to 10 
else 

write ( 6, *)' There are ',1,' elements.' 
write(4,*)l 
12=1 
endif 

elseif (string (4 : 6) .eq. ' 752' ) then 
ltmp=0 
nn (1) =0 

read(l, ' (A) ' ) string 
if (string (19 : 20) . eq. ' O') then 
ltmp=ltmp+l 

read (string (57 : 60) , ' (14) ' ) numb 
nn (ltmp+1) =numb+nn (ltmp) 
elseif (string(5: 6) .eq. ' -1' ) then 

write (6, *)' There are ' , ltmp, ' surfaces.' 
go to 20 

elseif (string (1 : 1 ). eq. ' ') then 
do i=l,4 

if (string ( (20*i-3) : 20*i) . ne . ' ') then 

1=1+1 

read (string ( (20*i-3) :20*i) , ' (14) ' )tr (1) 
endif 
enddo 
endif 
go to 30 
endif 
go to 15 
close (1) 

.Data from 'noddy' and 'elno' stored in 

nl - table of nodes making up the corresponding element 
x,y,z - nodal coordinate table 
rewind 2 
rewind 3 
do i=l,ll 

read (2, * ) na, x (i) ,y(i) ,z(i) 
enddo 
do i-1,12 

read (3, *)nk,nl (i, 1) ,nl (i,2) ,nl (i, 3) ,nl (i, 4) ,mat (i) 
enddo 

.Close 'noddy' and 'elno' for good 
close (2) 
close (3) 
ncount=0 

write (6, *) ' Be patient #*!?/#0*!!!' 
do ne=l,12 

.Store edge info in an integer array 'tab' after checking for 
overlap. The subroutine 'compare' is the heart of the program, 
call compare (tab) 

.Commenting out the following statement causes a speedup of 250% 

write (6, *)' Processed element no. ',ne,' :edge count= ' , ncount 
enddo 



ik/users / +research/+arnd/+ideas /proc 


PAGE 3 


07/26/91 3:16 PM 


write (6, *)' Edge count = ' ,ncount 
rewind (4) 

write ( 6, * ) ' No . of dielectric layers?' 

read (5, * ) itemp 

read (4,*) (xtmp, i=l , itemp) 

read(4,*) (nj,i=l,2) 

do i=l , (6*12) 

read (4 , * ) nk, edgv (i, 1) , edgv (i, 2) , edgv (i, 3) , n j 
enddo 
close (4 ) 
close (7) 
nsurf=0 

if (ivar.eq.l) then 
do i=l , ltmp 

write (6, *)' Processing surface ',i,' for on-surface edges.' 
do j=nn (i) +1, nn (i+1) 
do k=nn (i) +l,nn (i+1) 
do l=l,ncount 

if (tab (1, 1) . ne . 0) then 

if ( (tr ( j ) .eq. tab (1, 1) ) . and . (tr (k) . eq. tab (1, 2) ) ) then 
write (8, *) 1 
tab (1, 1) =0 
nsurf =nsurf+l 
go to 40 
endif 
endif 
enddo 
enddo 
enddo 
enddo 

write (6, *) 'Total no. of on-surface edges= ', nsurf 
write (6, *) 'Fem matrix to be solved is of order ', (ncount-nsurf ) 
else 
ncnt=0 

do ii=l,ltmp 
do ik=l,12 
nun (ik) =0 
enddo 

do i=nn (ii) +l,nn (ii+1) 
do j=l, 12 
do k=l , 4 

if (tr(i) .eq.nl (j, k) ) then 
nun ( j) =nun ( j) +1 
bc(j,nun(j) ) =nl (j,k) 
if (nun ( j) .eq. 3) then 

call edge ( j , be ( j , 1 ) , be ( j , 2 ) , be ( j , 3 ) ) 
ncnt=ncnt+l 
go to 80 
endif 
endif 
enddo 
enddo 
enddo 
enddo 

write (6, *)' There are ',ncnt,' surface elements.' 
endif 
close (8) 
stop 
end 
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compare' checks for all possible edges avoiding overlap and stores the 
info in an array 'tab'. The file 'edgy' contains the element no. 
and the six edge vectors corresponding to that element. 


subroutine compare (tab) 

integer nl (3000, 4) , mat (3000) 

real x (800) ,y (800) , z (800) 

integer nt (6, 2) , tab (4300, 2) 

common /bank/nl , x, y, z, ne, mat 

common /local/ncount 

nt (1, 1) =nl (ne, 1) 

nt (1, 2) =nl (ne, 2) 

nt (2, 1 ) =nl (ne, 1) 

nt (2, 2) =nl (ne, 3) 

nt (3, 1) =nl (ne, 1) 

nt (3, 2) =nl (ne,4) 

nt (4, 1) =nl (ne, 2) 

nt (4, 2) =nl (ne, 3) 

nt (5, 1 ) =nl (ne, 2) 

nt (5, 2) =nl (ne, 4) 

nt (6, 1) =nl (ne, 3) 

nt (6, 2) =nl (ne, 4) 

.Arranging node couplets in ascending order; a personal choice 
do i=l,6 

if (nt (i, 1) . gt . nt (i, 2) ) then 
ntmp=nt ( i , 1 ) 
nt (i, 1) =nt (i, 2) 
nt (i, 2) *ntmp 
endif 
enddo 

.A brute force search for overlapping edges 
l=ncount 
do ii=l , 6 

if (ne.eq.l) go to 32 
do i=l,ncount 

if ( (tab (i, 1) .eq.nt (ii, 1) ) .and. (tab (i, 2) .eq.nt (ii, 2) ) ) then 
write (4, *) ne, i, nt (ii, 1) , nt (ii, 2) , mat (ne) 
go to 35 
endif 
enddo 
1 = 1+1 

tab (1, 1) =nt (ii, 1 ) 
tab (1, 2) =nt (ii, 2) 

write(4,*)ne,l,nt(ii,l),nt(ii,2) ,mat (ne) 

write (7, *) l,x (nt (ii, 2) ) -x (nt (ii, 1) ) ,y (nt (ii, 2) ) -y (nt (ii, 1) ) , 

1 z (nt (ii, 2) ) -z (nt (ii, 1) ) 

enddo 
ncount=l 
return 
end 


determines edges corresponding to the surface element 

subroutine edge (m, jl, j2, j 3 ) 
integer edgv (18000, 3) , tmp (3) 
common /dbase/edgv 
do i j= (6*m) -5, (6*m) 

if ( jl . eq, edgv (i j , 2) ) then 
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if ( j2 . eq. edgv (i j, 3) ) then 
trap ( 1 ) =edgv ( i j ; 1 ) 
elseif ( j3 .eq. edgv (i j, 3) ) then 
tmp { 2 ) =edgv ( i j , 1 ) 
endif 

elseif ( j2.eq.edgv(i j,2) ) then 
if ( j3 . eq. edgv (i j , 3) ) then 
tmp (3) =edgv (i j, 1) 
endif 
endif 
enddo 

write (8,*)m,jl,j2, j3, tmp 

return 

end 
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.Computes the eigenvalues of an inhomogeneous cavity having arbitrary 
shape using edge elements. 

program fern 

integer edna (15000 ) , gnn (15000,2) , tab (6, 3) , sed (2000 ) , edst (3000,2) 

1 , mat (2500 ) 

character bit(6)*l 

real xyz (3000,3) ,al (6, 6) , x (500) , y (500) , z (500) , bl (6, 6) , eps (10) 
double precision a (700, 700) , b (700, 700) , alf r (700) , alfi (700) 

1 , beta (700) 

logical matz 

common /bank/edna, gnn, xyz, eps, mat 
common /mess/x,y,z 
.Reading in info 

write (6, *)' Number of edges' 
read (5, *) ned 

write (6, *) 'Number of on-surface edges' 
read(5,*)nes 

write (6, * ) ' If inhomogeneous, enter 1' 
read (5, *) ivar 
if (ivar.eq.l) then 

write (6, *)' Number of different permittivities' 
read (5, * ) nl 
endif 

open (1, file=' eglob' ) 
open (2, file=' edgy' ) 
open (3, file=' enoddy' ) 
open (4, file=' eigl' ) 
open (7, f ile=' esurfed' ) 
read(l,*) (eps (i) , i=l, nl) 
read (1, *) nn 
do i=l,nn 

read (3, * ) k,x(i),y(i),z(i) 
enddo 

write (6, *)' Finished reading in nodes' 
read(l, *) nel 
do i=l,6*nel 
itemp=l+ (i/6 . ) 

read (1, *) elm, edna (i) , gnn (i, 1) , gnn (i, 2) , mat (itemp) 
enddo 

write (6, *) eps (mat (1 ) ) , eps (mat (47) ) , eps (mat (48) ) , eps (mat (74 ) ) 
write (6, *) 'Finished reading in elements' 
do i=l,ned 

read (2, *) k,xyz (i, 1) ,xyz (i,2) ,xyz (i, 3) 
enddo 

write (6, *)' Finished reading in edges' 

read(7,*) (sed (i) , i=l, nes) 

close (1) 

close (2 ) 

close (3) 

close (7) 

.The following program segment generates the FEM matrix and imposes 
the boundary condition for the on-surface nodes simultaneously. 

It compares the on-surface edges stored in the array ' sed' and 
stores the complement in 'a' and 'b'. It also keeps a record of 
the new edge numbers and their corresponding old edge numbers in 
the pointer array 'edst'. 
net=ned-nes 
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do i=l,net 
do j=l,net 
a (i, j ) =0 . dO 
b (i, j ) =0 .dO 
enddo 
enddo 

write ( 6, *)' Generating FEM matrix' 
do i=l ,3000 
edst (i, 1) =0 
edst (i, 2) =0 
enddo 
nptrx=l 
do i=l,nel 

.Generate the matrix elements for the corresponding finite element 
and store the result in 'al' and 'bl' . 'tab' is a pointer array 
which stores the node combinations and local and global edge nos. 
call crux (i, al, bl , tab) 
do ij=l, 6 
bit (i j) =' 0' 
enddo 

.Run a check for the on-surface edges in the corresponding finite 
element. The character array 'bit' sets a flag for each on-surface 
edge in the finite element, 
do 100 j=l, 6 
do ichk=l,nes 

.If the element edge is an on-surface edge, set bit (i) ='l' 
if (tab ( j, 3) . eq. sed (ichk) ) then 
bit(j)='l' 
go to 100 
endif 
enddo 
enddo 

.If bit(i)= '1', forget it. If bit(i) ='0', check to see whether 
the edge has occurred previously; if not, increment 'nptrx' and 
store the new edge no. and corresponding old edge no. in 'edst' 
do j=l,6 

if (bit ( j) .eq. ' 0' ) then 

if (edst (tab ( j , 3) , 2) . ne . tab ( j , 3) ) then 
edst (tab ( j, 3) , 1 ) =nptrx 
edst (tab ( j , 3) , 2) =tab ( j, 3) 
nptrx=nptrx+l 
endif 
endif 
enddo 

.Only the edges inside the body are stored in the final array, 

' a' and ' b' . 
do j=l,6 
do k=l , 6 

if ( (bit ( j ) . eq. ' 0' ) . and. (bit (k) . eq. ' 0' ) ) then 
a (edst (tab ( j, 3) , 1) , edst (tab (k, 3) , 1) ) = 

1 a (edst (tab ( j , 3) , 1 ) , edst (tab (k, 3) ,1) ) +dble (al ( j , k) ) 
b (edst (tab ( j, 3) , 1) ,edst (tab (k, 3) , 1 ) ) = 

1 b (edst (tab ( j, 3) , 1) , edst (tab (k, 3) , 1) ) +dble (bl ( j, k) ) 
endif 
enddo 
enddo 
enddo 

write (6, *)' Finished FEM matrix generation' 
write (6, *) 'Matrix is of order ',nptrx-l 
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if (net . ne . (nptrx-1) ) then 

write (6, *)' Error in matrix assembly' 
stop 
endif 
ir=0 

do i=l, nptrx-1 
do j=l, nptrx-1 

if ( (b(i, j) -b(j,i) ) .gt. .000001) then 
write (6,*)b(i, j) ,b(j,i) ,i, j 
ir=l 
endif 
enddo 
enddo 

if (ir.eq.O) then 

write (6, *)' Symmetry test passed' 
else 

write (6, *)' Symmetry test failed' 
endif 

Solving the generalised eigenvalue problem: Ax=lambda*Bx 
epsl=. 000001 
matz=. false. 

call qzhes (700, net , a, b,matz, z) 
write (6, *)' Step 1 done' 
call qzit (700 , net , a, b, epsl ,matz, z, ier) 
write (6, *)' Step 2 done' 
write (6, * ) ' error= ',ier 
call qzval (700, net , a, b, alf r, alf i, bet a, mat z, z) 
write (6, *) ' Step 3 done' 
do i=l, nptrx-1 

alfr(i)=alfr(i) /beta (i) 
enddo 

A bubble sort through the eigenvalue spectrum 
do i=l, net-1 

do j=net,i+l,-l 

if (alf r ( j ) . It . alf r ( j — 1 ) ) then 
beta ( j ) =alf r ( j ) 
alfr( j)=alfr( j-1) 
alf r ( j-1 ) =beta ( j ) 
endif 
enddo 

if (alf r (i) . ge . 0 . ) then 
write (4, *) sqrt (alfr (i) ) 
endif 
enddo 
close (4 ) 
stop 
end 

subroutine crux (1, a, b, tab) 

integer edna (15000) , gnn (15000,2) , tab (6, 3) ,mat (2500) 

real xyz(3000,3),x(500),y(500),z(500),a(6,6),b(6,6),f(6,3), 

1 g (6, 3) , tmp (3) , eps (10) 

common /bank/edna, gnn, xyz, eps, mat 
common /mess/x,y,z 

common /local/sumx, sumy, sumz, xx, yy, zz, xy, yz, zx 
lv=6* (1-1) 
do j=l,6 

tab ( j, 1) =gnn (lv+ j, 1) 

tab ( j, 2) =gnn (lv+ j, 2) \ 
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tab ( j, 3) =edna (lv+ j) 
enddo 

.Sorting the array 'tab' arranges it according to local node nos. so 
that the array looks like the one in file 'input', 
call sort (tab) 

.'calc' calculates some data corresponding to the element 
call calc (tab (1, 1 ) , tab (1, 2) , tab (2,2) ,tab(3,2) ) 

.'volume' computes six times the volume of the tetrahedral element 
call volume (tab (1, 3) , tab (2,3) , tab (3, 3) , vol) 
do j=l,6 

call cross (x (tab (7- j , 1) ) , y (tab (7- j, 1) ) , z (tab (7- j , 1) ) , x (tab ( 

1 7- j, 2) ) , y (tab (7- j , 2 ) ) , z (tab (7- j, 2) ) , tmp) 

.'bj' stores the length of the 'j' th edge. 

b j=sqrt ( (xyz (tab ( j, 3) , 1) **2) + (xyz (tab (j,3),2)**2)+ (xyz 
1 (tab ( j, 3) , 3) **2) ) 

.The constant vectors, 'f' and 'g', are computed and the result 
stored in 'f' and 'g', respectively, 
do k=l , 3 

f ( j, k) =b j*tmp (k) /vol 
g ( j, k) =b j*xyz (tab (7- j, 3) ,k) /vol 
if (j.eq.2) then 
g(j/k)=-l.*g( j,k) 
elseif (j.eq.5) then 
f ( j,k)=-l.*f ( j,k) 
g( j,k)=-l.*g( j,k) 
endif 
enddo 
enddo 
vol=vol/6 . 
do j=l,6 
do k=l,6 

.'a' contains the elemnts of the 'a' matrix in the main program. 

This equals the volume integral over curl W_i . curl W_j and 
is the same as 4.*g_i*g_j*v 

a ( j, k) =4 . *dot ( j, k, g) *vol 

.The function 'fl' computes the elements of the 'b' matrix in the 
main program. This is taken directly from the formulation and 
equals the volume integral over W_i . W_j 
b( j, k) =fl ( j, k, f ,g) *eps (mat (1) ) *vol 
enddo 
enddo 
return 
end 

subroutine sort (tab) 
integer tab (6, 3) 
nc=l 

do ik=4,2,-l 

if (nc.eq.l) then 
i j=6 
else 
i j=ik 
endif 

do ii=l,ij-l 

do j=ij,ii+l,-l 

if (tab ( j , nc) . It . tab ( j-1, nc) ) then 
do jk=l,3 

call exchg (tab ( j , jk) , tab (j-1, jk) ) 
enddo 
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endif 

enddo 

enddo 

nc=2 

enddo 

if (tab (5, 2) . It . tab (4, 2) ) then 
do jk=l,3 

call exchg (tab (5, jk) , tab (4 , jk) ) 
enddo 
endif 

call exchg (tab (5, 1) , tab (5, 2 ) ) 

return 

end 

subroutine volume (jl, j2, j3,v) 
integer edna (15000) ,gnn (15000,2) ,mat (2500) 
real xyz (3000, 3) , bl (3) , b2 (3) , b3 (3) ,eps (10) 
common /bank/edna, gnn, xyz, eps,mat 
do it=l,3 

bl (it) =xyz ( jl, it) 
b2 (it) =xyz ( j2, it) 
b3 (it) =xyz ( j3, it) 
enddo 

v=abs (bl (1) * ( (b2 (2) *b3 (3) ) - (b2 (3) *b3 (2) ) ) -bl (2) * ( (b2 (1) *b3 
1 (3) ) - (b2 (3) *b3 (1) ) ) +bl (3) * ( (b2 (1) *b3 (2) ) - (b2 (2) *b3 (1) ) ) ) 

return 
end 


subroutine cross (xl,yl, zl,x2,y2, z2, tempo) 

real tempo (3) ’ 

tempo (1) =yl*z2-y2*zl 

tempo (2) =zl*x2-z2*xl 

tempo (3) =xl*y2-x2*yl 

return 

end 


subroutine calc ( jl, j2, j3, j4) 
real x (500) ,y (500) , z (500) 
common /mess/x,y,z 

common / local /sumx, sumy, sumz , xx, yy, zz, xy, yz, zx 
sumx=x ( jl ) +x ( j 2 ) +x ( j 3 ) +x ( j4 ) 
sumy=y ( j 1 ) +y ( j 2 ) +y ( j 3 ) +y ( j 4 ) 
sumz=z (jl)+z(j2)+z(j3)+z(j4) 

xx=sumx*sumx+x (jl)*x(jl)+x(j2)*x(j2)+x(j3)*x(j3)+x(j4)*x(j4) 
yy=sumy*sumy+y (jl)*y(jl)+y(j2)*y(j2)+y(j3)*y(j3)+y(j4)*y(j4) 
zz=sumz*sumz+z (jl)*z(jl)+z(j2)*z(j2)+z(j3)*z(j3)+z(j4)*z(j4) 
xy=sumx*sumy+x (jl)*y(jl)+x(j2)*y(j2)+x(j3)*y(j3)+x(j4)*y(j4) 
yz=sumy*sumz+y (jl)*z(jl)+y(j2)*z(j2)+y(j3)*z(j3)+y(j4)*z(j4) 
zx=sumz*sumx+z (jl)*x(jl)+z(j2)*x(j2)+z(j3)*x(j3)+z(j4)*x(j4) 
return 
end 

real function dot ( jl, j2,vec) 
real vec(6,3) 

dot= (vec ( jl , 1 ) *vec ( j2, 1 ) ) + (vec ( jl, 2) *vec ( j2, 2) ) + (vec(jl,3) *vec 
(j2,3) ) 
return 
end 
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real function fl ( jl, j2, f ,g) 

real f(6,3),g(6,3), tmpl (3) , tmp2 (3) 

common /local/sumx, sumy, sumz,xx,yy, zz,xy,yz, zx 

call cross (f(jl,l),f(jl,2),f(jl,3),g(j2,l),g(j2,2),g(j2,3), tmpl) 
call cross (f(j2,l),f(j2,2),f(j2,3),g(jl,l),g(jl,2),g(jl,3), tmp2) 
terml=dot ( jl, j2, f ) 

term2= ( (tmpl (1 ) +tmp2 (1 ) ) *sumx+ (tmpl (2) +tmp2 (2 ) ) *sumy+ 

1 (tmpl (3) +tmp2 (3) ) *sumz) /4 . 

term3= (g(jl,2)*g(j2,2)+g(jl,3)*g(j2,3)) *xx+ (g(jl,l)*g(j2,l)+ 

1 g( jl,3)*g( j2,3) ) *yy+ (g ( jl, 1) *g( j2, 1) +g( jl,2) *g( j2, 2) ) *zz 

2 - (g ( jl, 1) *g ( j2, 2) +g ( jl, 2) *g ( j2, 1) ) *xy- (g ( jl, 1) *g ( j2, 3) + 

3 g( jl, 3) *g( j2, 1) ) *zx- (g(jl,2)*g(j2,3)+g(jl,3)*g(j2,2)) *yz 
f l=terml+term2+ (term3/20 . ) 

return 

end 


subroutine exchg(jl,j2) 

ntmp= jl 

jl = j2 

j2=ntmp 

return 

end 
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rocesses a universal file obtained from IDEAS and stores the nodal info. 
Ihecks the surface nodes and imposes suitable boundary conditions on them. 

itores the node numbers and respective nodal coordinates in 'noddy' 
itores the element numbers and corresponding nodes in 'elno' 
itorage limit: 4000 nodes, 25000 elements 

program unv_f ile_processor 

character string*80, yastrn*40, unv*20 

integer nl (25000,4) ,nn(50) , nun (25000) , be (25000, 3) 

integer count, tr (25000) , tmp 

real x ( 4000) , y (4000) , z (4000 ) 

write (6, *) ' 1) Cavity problem 2) Scattering problem' 
read (5, *) ivar 

write (6, *) 'Name of universal file ?' 

read (5, ' (A) ' ) unv 

open (1, file=unv) 

open (2, file=' noddy' ) 

open (3, file=' elno' ) 

open (8, f ile=' surfed' ) 

Processing universal file for info on various parameters 
read (1, ' (A) ' ) string 
1=0 

if (string (4 : 6) .eq. ' 151' ) then 
write (6, *)' Encountered header' 
go to 15 

Processing nodal info; data stored in 'noddy' 
elseif (string (4 : 6) .eq. ' 15') then 
read(l, ' (A) ' ) string 
if (string (5 : 6) . ne . ' -1' ) then 

write (2, *) string(7 : 10) , ' ', string (41 : 53) , ' ', string (54 : 66) , ' ' 

1 , string(67 : 79) 

1=1+1 
go to 5 
endif 
11=1 

write (6, *)' There are ',1,' nodes.' 
write (8, *) 1 

Processing element info; data stored in 'elno' 
elseif (string (4 : 6) . eq. ' 71') then 
read(l, ' (A) ' ) string 
if (string(5: 6) .ne. ' -1' ) then 
read(l, ' (A) ' )yastrn 

write (3, *) string (7 : 10) , ' ' ,yastrn (7 : 10) , ' ' , yastrn (17 : 20) , 

1 ' ', yastrn (27 : 30) , ' ' , yastrn (37 : 4 0) 

1=1+1 
go to 10 
else 

write (6, *) 'There are ',1,' elements.' 
write (8, *) 1 
12=1 
endif 

elseif (string (4 : 6) .eq. ' 752' ) then 
ltmp=0 
nn (1 ) =0 

readd, ' (A) ' ) string 
if (string(19:20) .eq. ' O') then 
ltmp=ltmp+l 
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read (string (57:60),' (14)') numb 
nn (ltmp+1 ) =numb+nn (ltmp) 
elseif (string (5 : 6) . eq. ' -1' ) then 

write (6, *) 'There are ' , ltmp, ' surface divisions.' 
go to 20 

elseif (string (1 : 1 ). eq. ' ') then 

do i=l,4 

if (string ( (20*i-3) : 20*i ). ne . ' ') then 

1=1+1 

read (string ( (20*i-3) :20*i) , ' (14) ' ) tr (1) 
endif 
enddo 
endif 
go to 30 
endif 
go to 15 
close (1 ) 

.Data from 'noddy' and 'elno' stored in 

nl,..,n4 - table of nodes making up the corresponding element 

x,y,z - nodal coordinate table 

rewind 2 
rewind 3 
do i=l ,11 

read (2, *)na,x(i),y(i),z(i) 
enddo 
do i=l , 12 

read (3, * ) nk, nl (i, 1 ) , nl (i, 2) , nl (i, 3) , nl (i, 4) 
enddo 

.Close 'noddy' and 'elno' for good 
close (2) 
ncnt=0 
do i=l,12 
do j=l,3 
be (i, j ) =0 
enddo 
enddo 

if (ivar.eq.l) then 
do i=l,ltmp 

write (6, *)' Processing surface ',i,' for on-surface nodes.' 
do j=nn (i) +2, nn (i+1) 

if (x (tr (nn (i) +1) ) .eq.x (tr ( j) ) ) then 
f lagx=0 - 
else 

f lagx=l . 
endif 

1 if (y (tr (nn (i) +1) > .eq.y (tr ( j) ) ) then 

f lagy=0 . 
else 

f lagy=l . 
endif 

if (z (tr (nn (i) +1) ) .eq. z (tr ( j) ) ) then 
f lagz=0 . 
else 

f lagz=l . 
endif 
enddo 

do k=nn (i) +1, nn (i+1) 
if (flagx.eq. 0 . ) then 
be (k, 1) =1 
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elseif (f lagy . eq. 0 . ) then 
be (k, 2 ) =1 

elseif ( f lagz . eq. 0 . ) then 
be (k, 3) =1 
endif 
enddo 
enddo 

do i=l/ 1-1 

do j=l,i+l,-l 

if (tr (j) .lt.tr (j-1) ) then 
call exchg (tr ( j ) , tr ( j-1 ) ) 
do k=l , 3 

call exchg (be ( j, k) , be ( j-1, k) ) 
enddo 
endif 
enddo 
enddo 
nuk=0 . 
do i=l , 1 

if (i.eq.l) then 
tmp=0 
else 

tmp=tr (i-1) 
endif 

if (tr (i) .ne.tmp) then 
count=l 

if (tr (i+count) . eq. tr (i) ) then 
count =count+l 
go to 50 
endif 
ntl=l 
nt2=l 
nt3=l 

do j=l, count 

ntl=ntl*bc (i+ j-1, 1) 
nt2=nt2*bc (i+ j-1, 2) 
nt3=nt3*bc (i+ j-1, 3) 
enddo 

write (8 , * ) tr (i) , ntl, nt2, nt3 
ncnt=ncnt+l 

nuk=nuk+ (1-ntl) + (l-nt2) + (l-nt3) 
else 

go to 100 
endif 
enddo 

write (6, *) 'No. of on-surface nodes = ' , nent 
write (6, *) 'No. of unknowns = ',3*ll-nuk 
else 

do i=l,ltmp 

write (8, * ) i, nn (i+1 ) 
enddo 

do ii=l,ltmp 
do ik=l,12 
nun (ik) =0 
enddo 

do i=nn (ii) +1, nn (ii+1) 
do j=l,12 
do k=l,4 

if (tr (i) .eq.nl (j, k) ) then 
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nun ( j ) =nun ( j ) +1 

be ( j, nun< j) ) =nl ( j, k) 

if (nun ( j ) . eq. 3) then 

write (8,*) j , be ( j , 1 ) , be ( j , 2 ) , be ( j , 3 ) 
ncnt=ncnt+l 
go to 200 
endif 
endif 
enddo 
enddo 
enddo 
enddo 

write (6, *)' There are ',ncnt,' surface elements.' 
endif 
close (3) 
close (8) 
stop 
end 


subroutine exchg(jl,j2) 

ntmp= jl 

j 1 = j 2 

j2=ntmp 

return 

end 
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This program computes the eigenvalues of a rectangular cavity 
using divergenceless nodal basis and tetrahedral elements. 

Storage limit: 500 nodes, 2000 elements, 300 surface nodes 
1500 matrix size 
To run program on Apollos : 

Need to run 'procn.obj' first; 'noddy', 'elno' and 'surfed' are created, 
ftn ' filename' 

bind ' filename' . bin /progs/math/naas/eispack . lib -b ' filename' . obj 
' filename' . obj 

program fem_nodal_basis_3d 
logical matz 

integer elm (2000 , 4 ) , surf (300) ,bc(300,3) ,tab(4) 
real al(12,4),bl(12,4),x(500),y(500),z(500) 

double precision a (1500,1500) , b (1500, 1500) ,alfr (1500) , alfi (1500) 

1 ,beta(1500) 

common /dbase/x, y, z, elm 
common /var/s, vol 

write (6, *) 'No. of on-surface nodes' 
read ( 5, * ) nsn 

write (6, *) 'Penalty factor' 

read(5, *) s 

open ( 1 , f ile=' noddy' ) 

open (2, file=' elno' ) 

open (3, file=' surfed' ) 

open (4 , f ile=' eig' ) 

read(3,*)nn 

nn2=2*nn 

read (3, *) ne 

do i=l,nn 

read(l, *)nk,x(i) ,y(i) ,z(i) 
enddo 
do i=l,ne 

read (2, * ) nk, elm (i, 1 ) , elm (i, 2) , elm(i, 3) , elm(i, 4) 
enddo 

do i=l,nsn 

read (3, * ) surf (i) , be (i, 1 ) , be (i, 2) , be (i, 3) 
enddo 

write (6, *)' Finished reading in data' 
close (1 ) 
close (2) 
close (3) 


do iv=l , 3 
do i=l,ne 

call crux (i, iv, al , bl, tab) 
do j=l , 12 
do k=l , 4 

nt= (iv-1) *nn 
if (j.le.4) then 

a (tab ( j ) , (nt+tab (k) ) ) =a (tab ( j ) , (nt+tab (k) ) ) +dble (al ( j, k) ) 
b (tab ( j) , (nt+tab (k) ) ) =b (tab ( j) , (nt+tab (k) ) ) +dble (bl ( j, k) ) 
elseif ( ( j . gt . 4 ) . and. ( j . le . 8) ) then 
a ( (nn+tab ( j-4 ) ) , (nt+tab (k) ) ) =a ( (nn+tab ( j — 4 ) ) , (nt+tab (k) ) ) 
1 +dble (al ( j, k) ) 

b ( (nn+tab ( j-4 ) ) , (nt+tab (k) ) ) =b ( (nn+tab ( j-4 ) ) , (nt+tab (k) ) ) 
+dble (bl ( j , k) ) 
else 
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a ( (nn2+tab ( j-8 ) ) , (nt+tab (k) ) ) =a ( (nn2+tab ( j-8) ) , (nt+tab (k) ) ) 
1 +dble (al < j, k) ) 

b ( (nn2+tab (j-8) ) , (nt+tab (k) ) ) =b ( (nn2+tab ( j-8 ) ) , (nt+tab (k) ) ) 
1 +dble(bl ( j,k) ) 

endif 
enddo 
enddo 
enddo 
enddo 


write (6, *) 'Applying be' 
do j=l,nsn 

if (be ( j, 1) .eq. 0) then 
a (surf (j) , surf (j) )=-l.d0 
endif 

if (be ( j, 2) .eq. 0) then 

a (nn+surf ( j ) , nn+surf ( j ) ) =-l . dO 
endif 

if (be ( j, 3) .eq. 0) then 

a (nn2+surf ( j ) , nn2+surf ( j ) ) =-l . dO 
endif 
enddo 

write ( 6 ,*)' shifting rows' 

nptr=3*nn+l 

do j=l, (3*nn) 

if (f loat ( j/15) . eq. ( j/15 . ) ) then 
write (6, *) j, nptr 
endif 

if (a ( j, j) .eq. -1 .dO) then 
nptr=nptr-l 

if (j.gt.nptr) go to 30 
if (a (nptr, nptr) . ne . -1 . dO) then 
do k=l,nptr 

call exchg (a ( j , k) , a (nptr , k) ) 
call exchg (b ( j , k) , b (nptr, k) ) 
enddo 

do k=l,nptr 

call exchg (a (k, j) , a (k, nptr) ) 
call exchg (b (k, j ), b (k, nptr) ) 
enddo 
else 

go to 20 
endif 
endif 
enddo 

ir=0 

write (6, *) 'Matrix generation complete; matrix of order ',nptr 
do i=l,nptr 
do j=l,nptr 

if ( (a (i, j ) -a ( j , i) ) .gt . . 000001 ) then 
ir=l 

write (6, *) a (i, j),a(j,i),j,i 
endif 
enddo 
enddo 

if (ir.eq.l) then 

write (6, *)' Symmetry test failed for A' 
else 
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write (6, *)' Symmetry test passed for A' 
endif 


epsl= . 00001 
matz= . false . 

call qzhes (1500, nptr, a, b,matz, z) 
write (6, *)' Step 1 done' 

call qz it (1500, nptr, a, b, epsl , matz, z, ier) 
write (6, *)' Step 2 done' 
write (6, *) ' error= ' , ier 
if (ier.ne.O) then 

write ( 6, *)' Erroneous computation; check FE matrices' 
stop 
endif 

call qzval (1500, nptr, a, b, alf r, alf i , beta, mat z , z) 
write (6, *)' Step 3 done' 
do i=l,nptr 

alfr(i)=alfr(i) /beta (i) 
enddo 

do i=l,nptr-l 

do j=nptr, i+1 , -1 

if (alf r ( j ) . It . alf r ( j -1 ) ) then 
beta ( j ) =alf r ( j ) 
alfr ( j)=alfr ( j-1) 
alf r ( j-1 ) =beta ( j ) 
endif 
enddo 

if (alfr (i) .gt . 0. ) then 
write (4, *) sqrt (alfr (i) ) 
endif 
enddo 
close (4 ) 
stop 
end 


subroutine crux (1, iv, a, b, tab) 

integer elm (2000, 4 ), tab (4 ) 

real a(12,4) , b (12, 4 ) ,x(500) ,y(500) , z(500) 

common /dbase/x, y, z, elm 

common /var/s,vol 

do j=l,4 

tab ( j) =elm(l, j ) 
enddo 

call calc (iv, tab (1 ) , tab (2) , tab (3) , tab (4 ) ) 
call volume (tab (1 ) , tab (2) , tab (3) , tab (4 ) , vol) 
do j=l,12 
do k=l,4 

call matx (iv, j , k, a ( j, k) , b ( j , k) ) 
enddo 
enddo 
return 
end 

subroutine volume (jl, j2, j3, j4,v) 
integer elm(2000,4) 
real x (500) ,y (500) , z (500) 
common /dbase/x, y, z, elm 

Vl=(x(jl)-x(j4))*((y(j2)-y(j4))*(z(j3)-z(j4))-(y(j3)-y(j4))*(z(j 
2) -z ( j4) ) ) 
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V2=(y(j4)-y(jl))*((x(j2)-x(j4))*(z(j3)-z(j4>)-(x(j3)-x(j4))*(z(j 
1 2 ) -z ( j4 ) ) ) 

v3= ( z ( jl) -z ( j4 ) ) * ( (x ( j2) -x ( j4 ) ) * (y ( j3) -y ( j4 ) ) - (x ( j3) -x ( j4 ) ) * (y ( j 

1 2) -y ( j 4 ) ) ) 

v=abs (vl+v2+v3) /6 . 

return 

end 

subroutine calc (iv, jl, j2, j3, j 4 ) 
integer elm(2000,4) 

real x (500 ) , y ( 500 ) , z (500 ) , p (4 ) , q (4 ) , r ( 4 ) 
common /dbase/x, y, z, elm 
common /local/p, q,r 

P(l) = (y(j2)-y(j4))*(z(j3)-z(j4))-(y(j3)-y(j4))*(z(j2)-z(j4)) 

P(2) = (y(j3)-y(j4) )*(z(jl)-z(j4) ) — (y(jl) — y ( j4 ) )*(z(j3)-z(j4) ) 
p(3) = (y(jl)-y(j4) )*(z(j2)-z(j4) )-(y(j2)-y(j4) )* (z ( jl ) — z ( j4 ) ) 

P(4) = (y(j3)-y(jl))*(z(j2)-z(jl))-(y(j2)-y(jl))*(z(j3)-z(jl)) 
q(l) = (z (j2) -z (j4) ) * (x(j3) -x(j4) ) - (z ( j 3 ) -z (j4) ) * (x(j2) -x(j4) ) 
q(2) = (z(j3) -z ( j 4 ) ) * (x ( jl) -x ( j4) )-(z(jl)-z(j4) )* (x(j3)-x(j4) ) 
q(3) = (z (jl) -z (j4) ) * (x(j2) -x(j4) ) -(z { j 2 ) -z { j 4 ) ) * (x ( jl ) -x { j4 ) ) 
q(4) = (z (j3) -z (jl) ) * (x(j2) -x(jl) ) -(z (j2) -z (jl) ) * (x(j3) -x ( jl) ) 
r(l) = (y(j3) -y(j4) )* (x(j2)-x(j4) ) — (y {j2) — y ( j 4 ) )* (x(j3)-x(j4) ) 
r (2 ) = (y ( jl ) -y ( j4 ) ) * (x ( j3) -x ( j4 ) ) - (y ( j3) -y ( j4 ) ) *(x ( jl ) -x ( j4 ) ) 

r (3 ) = (y ( j2 ) -y ( j4 ) ) * (x ( jl ) -x < j4 ) ) - (y ( jl ) -y (j4) )* (x(j2)-x(j4) ) 
r(4) = (y(j2)-y(jl) ) * (x ( j3) -x ( jl ) ) - (y { j3) -y ( jl) )* (x(j2)-x( jl) ) 
Checking routine follows; can be deleted if the user feels confident 
enough of the matrix generation process 
tl=0 . 
t2=0 . 
t3=0 . 

do ijk=l,4 
tl=tl+p (i jk) 
t2=t2+q(i jk) 
t3=t3+r (i jk) 
enddo 

if (tl.gt. .000001) then 
write (6, *)tl,t2,t3 
endif 
return 
end 


subroutine matx (m, j , k, f a, fb) 
real p(4) ,q(4) , r (4) 
common /local/p, q,r 
common /var/s,vol 
tt=0 . 

if (m.eq.l) then 
if (j.le.4) then 

t=(s*p( j)*p(k) ) + (q( j) *q(k) ) + (r( j) *r <k) ) 
if (j.ne.k) then 
tt=. 05 
else 
tt=. 1 
endif 

elseif ( ( j . gt . 4 ) . and. ( j . le . 8) ) then 
j j= j-4 

t= (s*q( j j) *p (k) ) - (p( j j) *q(k) ) 
else 
j j=j-8 
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t=(s*r(jj)*p(k))-lp(jj)*r(k)) 

endif 

elseif (m.eq.2) then 
if (j.le.4) then 

t=(s*p(j)*q(k) )-(q(j)*p(k) ) 
elseif ( ( j . gt . 4) . and. ( j . le . 8) ) then 

j j=j_4 

t=(p( j j) *P(M ) + (s*q( j j) *q(k) ) + (r ( j j) *r (k) ) 
if (jj.ne.k) then 
tt= . 05 
else 
tt=.l 
endif 
else 
jj=j-8 

t=(s*r(jj) *q(k) ) — ( q ( j j ) *r(k) ) 
endif 
else 

if (j.le.4) then 

t=(s*p(j) *r(k) ) - (r ( j) *p(k) ) 
elseif ( ( j . gt . 4 ) . and. ( j . le . 8) ) then 
j j=j-4 

t= (s*q ( j j ) *r (k) ) - (r ( j j) *q(k) ) 
else 
jj=j-8 

t=(p( j j) *p (k) ) + <q( j j) *q(k) ) + (s*r ( j j) *r (k) ) 
if (jj.ne.k) then 
tt=. 05 
else 
tt=.l 
endif 
endif 
endif 

f a=t/ (36 . *vol) 

fb=tt*vol 

return 

end 


subroutine exchg(xl,x2) 

tmp=xl 

xl=x2 

x2=tmp 

return 

end 
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